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1. Introduction. , 

This period included many scientific observations made with the Chandra Observatory. 

The results, as is well known, are spectacular. Fortunately, the HRMA performance continues 
to be essentially identical to that predicted from ground calibration data. The Telescope 
Scientist Team has improved the mirror model to provide a more accurate description to 
the Chandra observers and enable them to reduce the systematic errors and uncertainties 
their data reduction. We also have made considerable progress in improving the scattering 

model. 

There also has been progress hi the scientific program. At this time 58 distant clusters 
of galaxies have been observed. We are performing a systematic analysis of this rather large 
data set for the purpose of determining absolute distances utilizing the Sunyaev Zel dovich 
effect. These observations also have been used to study the evolution of the cluster baiyon 
mass function and the cosmological constraints which result from this evolution. 

2. Chandra Mirror Activities 

2.1. Monitoring and Program Support . , 

We continued to monitor the HRMA mirror performance in flight; during this period we 
did not find it necessary to update the basic mirror parameters or ray tracing models as a 
result of changes in flight performance. We contributed to NASA s Chandra A04 Research 
Announcement. We rewrote and updated the HRMA chapter of the ‘‘Chandra Proposers 
Observatory Guide”, which was sent to all the Chandra proposers with the NASA s Chandra 
AOs We also worked with CXC optical group for the ray-trace analysis, with the Chandra 
user support group to provide information needed for the Chandra users, and with the CXC 
data center to provide information for the CIAO database. 

2 2 Scattering 

' We continued our study of X-ray scattering from random rough surfaces, and established 
a new scattering model. 

The theory of scattering from random rough surfaces has been the subject of a number 
of books and many research papers. The modeling of X-ray scattering at grazing angles 
is difficult because of its short wavelength and the small angle between the wave and the 
surface. Almost all the current approaches use a common approximation the scattering 
angle is much smaller than the grazing incident angle. Most of the approaches also ignore 
the scattering asymmetry around the direction of specular reflection (scattering towards vs. 
away from the surface). Many of the present models also are limited to small perturbations 
of the ideal surface. These approximations are not adequate for many of the applications 
involving X-ray mirrors. 
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To fully evaluate the Chandra performance and analyze the scientific data, we need 
an accurate model of the X-ray scattering. Our new method models the X-ray scattering 
directly from first principles. An actual rough surface is (incompletely) described by its 
Power Spectral Density (PSD). For a given PSD, model surfaces with the same roughness as 
the actual surface are constructed by preserving the PSD amplitudes and assigning a random 
phase to each spectral component. Rays representing the incident wave are reflected from 
the model surface and projected onto a flat plane, which approximates the model surface, 
as outgoing rays and corrected for phase delays. The projected outgoing rays are then 
corrected for wave densities and redistributed onto an uniform grid where the model surface 
is constructed. The scattering is then calculated by taking the Fast Fourier Transform (FFT) 
of the resulting distribution. 

Our new method treats the reflection and scattering together, and consequently both 
depend upon the surface roughness. It does not require the approximation that the scattering 
angle is small compared to the grazing angle so that all the scattered rays can be traced 
accurately. This method is generally applicable. We have applied this method to the mirrors 
of the Chandra X-ray Observatory and have shown that the calculated scattering profile is 
as expected, including the Fraunhofer scattering patterns which result from the finite length 
of the surfaces. 

This work is still continuing. Next we will generate scattering tables, which are the 
tabulations of the scattering function S. Then we will use these scattering tables in our 
ray-trace model to simulate the CXO performance and compare it with the real results of 
the CXO, from both on-orbit observations and its ground calibrations. This method should 
be useful for other X-ray telescope missions as well. 

We have presented this work at the SPIE ’02 meeting. The paper describing this new 
method is published in the SPIE proceedings 4851-11 and is included as appendix A of this 
report. 

3. The Funtools Project 

The Funtools project arose out of conversations with astronomers about the decline in 
their software development efforts over the past decade. A stated reason for this decline is 
that it takes too much effort to master one of the existing FITS libraries simply in order to 
write a few analysis programs. This problem is exacerbated by the fact that astronomers 
typically develop new programs only occasionally, and the long interval between coding 
efforts often necessitates re-learning the FITS interfaces. 

As stated in our proposal for this past year, software development plans for the funtools 
software suite and the xpa messaging system covered the following areas: 

• Closer integration of funtools with the ds9 image display program. 

• Enhanced security of funtools and xpa code. 

• Development of firewall support in xpa. 
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• Provide elliptical and rectangular annuli support for funtools region masks. 


In all four areas we accomplished significant advances. 

1. A graphical parameter interface was developed in ds9 to support user-specified inputs 
to analysis routines. Macro expansion of funtools command lines also was enhanced so that 
various internal ds9 values could easily be passed to the funtools programs. As a result, 
funtools programs now are routinely executed from within ds9. These enhancements weie 
released in ds9 version 2.1 in April 2002. 

2. Code security was enhanced for both funtools and xpa. The insecure Unix sys- 
tem!) function was replaced by our secure launchQ routine. The latter explicitly utilizes 
fork () /exec () instead of ”sh -c” and thus is not susceptible to the usual shell-based attacks. 
The insecure Unix mktemp() routine was replaced by our secure mkrtempQ routine. The 
latter makes a unique temporary filename that is hard to guess and optionally opens it us- 
ing secure flags that prevent race conditions. These enhancements were released in funtools 

version 1.2 in April 2002. 

3 The xpa system now supports communication through firewalls. Using the new xpare- 
mote command, a local xpa-enabled program such as ds9 can register itself with a remote 
xpans name server set in xpa proxy mode. The local ds9 can then receive data and com- 
mands from behind a port-blocking firewall, since the socket connection is made securely 
from inside the firewall to the xpans on the outside. Network address translation (NA1) 
schemes also are handled properly. These enhancements were released in xpa version 2.1 in 


April 2002 

4. Elliptical and rectangular annuli regions were added to funtools/ds9 spatial masks. 
The panda (Pie AND Annulus) shape also was extended to ellipses and boxes. The syntax 


for these commands is shown below: , . . 

BOX xcenter ycenter xwl yhl xw2 yh2 ... xwn yhn (angle) 

BOX xcenter ycenter xwlo yhlo xwhi yhhi n= [number] (angle) 

ELLIPSE xcenter ycenter xwl yhl xw2 yh2 ... xwn yhn (angle) 

ELLIPSE xcenter ycenter xwlo yhlo xwhi yhhi n= [number] auto(angle) 

BPANDA xcen ycen angl ang2 nang xwlo yhlo xwhi yhhi nrad (ang) # box 

EPANDA xcen ycen angl ang2 nang xwlo yhlo xwhi yhhi nrad (ang) # ellipse 

These enhancements were first released in beta test versions of funtools starting in 

September 2002. 


For general information on xpa and funtools (including download information), please 

see: http:/ /hea- www.harvard.edu/RD/xpa/index.html 

http: / /hea- www .harvard.edu/RD/ funtools/ index, html 


4. Science Program. 

The Telescope Scientist guaranteed time is devoted to observations of distant galaxy 
clusters to study their evolution and to determine cosmological scale quantities using the 
Sunyaev Zel'dovich effect. Cooperative agreements with other scientists, including radio 
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observers, have been formed. Preliminary data processing has been performed upon the 58 
data sets received at this time. This processing includes correction for the changes in the 
ACIS gain, the identification and elimination of periods of high background, the removal of 
ACIS after-glow effects, the calculation of exposure maps and images, and the search for 
serendipitous sources. The superior angular resolution of Chandra shows more structure in 
galaxy cluster systems than we had expected based upon pre-Chandra observations; this 
significantly complicates the final data analysis. 

The different sensitivity of the radio and X-ray measurements to the outer parts of the 
cluster contributes to a basic uncertainty in the analysis of these data for the purpose of 
determining distances. We have begun a program of simulations of the X-ray and radio data 
to determine the sensitivity of the final results to the relatively poorly measured low surface 
brightness regions at large cluster radii. 

One very important result already has been obtained, and is in the process of final 
preparation for publication by Vikhlinin et al. We find that the cluster density for equivalent 
baryon masses strongly evolves between redshifts of 0.4-0.8 and the present epoch; the local 
density of mass clusters per unit co-moving volume is approximately 10 times that found at 
redshifts of 0.5. The density evolution is a strong function of the cosmological parameters ft 
and A. We find that the band in these parameters allowed by cluster evolution is comparable 
in scale to the bands allowed by SN la and the CMB, and that all three bands intersect in a 
small region near the currently commonly used values of ft = 0.3, A = 0.7. This may either 
be viewed as strong supporting evidence that there are not important systematic errors in 
the interpretation of the three types of data, or, if we accept the cosmological parameters 
as determined form the SN la and CMB data, then we have strong supporting evidence for 
models of cluster formation. 

Dr. Ping Zhao continued the study of stellar black holes. This includes the search for 
new black holes and the study of known black holes. Dr. Zhao is a principal investigator of 
a program to monitor the quiescent black hole X-ray novae (BHXN), and also is a principal 
investigator of a program to search for quiescent X-ray novae among the old optical novae. 

Dr. Zhao also is a co-investigator of the Chandra Multiwavelength Plane Survey Project 
(ChaMPlane) which is a project to identify a large sample of serendipitous X-ray sources in 
the galactic plane Chandra fields, in order to determine the populations of accretion-powered 
binaries in the Galaxy. He has supported these studies with optical observations using the 
MMT and Magellan telescopes. 
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A. Paper on X-ray Scattering Published in 2002 SPIE 4851-11 
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A new method to model X-ray scattering from 
random rough surfaces 

Ping Zhao and Leon P. Van Speybroeck 

Harvard-Smithsonian Center for Astrophysics 
60 Garden Street, Cambridge, MA 02138 U.S.A. 


ABSTRACT 

This paper presents a method for modeling the X-ray scattering from random rough surfaces. An actual rough 
surface is (incompletely) described by its Power Spectral Density (PSD). For a given PSD, model surfaces with 
the same roughness as the actual surface are constructed by preserving the PSD amplitudes and assigning a 
random phase to each spectral component. Rays representing the incident wave are reflected from the model 
surface and projected onto a flat plane, which approximates the model surface, as outgoing rays and corrected 
for phase delays. The projected outgoing rays are then corrected for wave densities and redistributed onto 
an uniform grid where the model surface is constructed. The scattering is then calculated by taking the Fast 
Fourier Transform (FFT) of the resulting distribution. This method is generally applicable and is not limited 
to small scattering angles. It provides the correct asymmetrical scattering profile for grazing incident radiation. 
We apply this method to the mirrors of the Chandra X-ray Observatory and show the results. We also expect 
this method to be useful for other X-ray telescope missions. 

Keywords: X-ray scattering, random rough surface, X-ray mirror, X-ray telescope, Chandra X-ray Observatory 


1. INTRODUCTION 

The study of scattering from random rough surfaces goes back at least to Rayleigh in 1887, 1 and has been 
investigated by many physicists and engineers. The problem has been the subject of numerous books, including 
the classic “The Scattering of Electromagnetic Waves From Rough Surfaces 1 ’ by Beckmann and Spizzichino 2 
and countless research papers. 3 Many people have also studied X-ray scattering at grazing angles as part of this 
subject. This problem is even more difficult because of the short wavelength (compared to the scale of the surface 
roughness) and the small angle between the wave propagation direction and the surface. Most approaches in 
the literature make the approximation that the scattering angle is much smaller than the incident grazing angle. 
Some of the treatments use the approximation that the surfaces are sufficiently “smooth” so that a low order 
expansion in the surface height errors is adequate, and consequently are limited in their applications. Many of 
the approaches can not obtain the scattering asymmetry around the direction of specular reflection (scattering 
towards vs. away from the surface). These approximations are not adequate for many of the applications 
involving X-ray mirrors. 

This new study of the century old problem is motivated by our direct involvement of the evaluation of 
the X-ray mirror performance aboard the Chandra X-ray Observatory (CXO) - the NASA’s third great space 
observatories now have been successfully operated for three years and have brought us fruitful scientific results 
with many exciting discoveries. A major achievement of the CXO compared to previous X-ray missions is its 
unprecedented spatial resolution (< 0.5" FWHM). This is mainly due to the design and manufacture of its 
X-ray mirrors. These mirrors are the largest, most precise grazing incidence optics ever built. At 0.84-m long 
and 0.6 - 1.2-m in diameters, the surface area of each mirror ranging from 1.6 to 3.2 square meters. They were 
polished to the highest quality ever achieved for any X-ray mirrors of this size. The surface roughness of these 
mirrors is comparable to or less than the X-ray wavelengths in the 0.1-10 keV band over most of the mirror 
surfaces. 
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However, the mirrors are not perfect, and consequently there are still small amount of scattered X-rays. We 
need an accurate model of the X-ray scattering to fully evaluate the CXO performance and analyze the scientific 
data. We have built a raytrace computer model to simulate the CXO performance. In the current raytrace 
model, the reflection and scattering are treated as separate effects. Each ray hit the mirror is reflected according 
the reflectivity of the surface geometry and the low frequency surface map, but not the high frequency roughness, 
therefore the effect of the surface roughness on the reflection efficiency is lost. The scattering is treated with 
small (scattering) angle approximation. 

Our new method treats the reflection and scattering together, and consequently both depend upon the 
surface roughness. It does not require the approximation that the scattering angle is small compared to t le 
grazing angle so that all the scattered rays can be traced accurately. 

2. POWER SPECTRAL DENSITY OF ROUGH SURFACES 

A rough surface is described, statistically, by its surface Power Spectral Density (PSD) as a function of the 
surface spatial frequency, f. Consider a 1-dimensional surface with length L and surface height (ie. deviation 
from a perfectly flat surface): 2 = h(x), which assumes the value s (-00 < 2 < 00). Its PSD is dehned as. 


PSD(f) = 2 Wi(f) = 7 


,L/2 

L ,2 


e ,2nxf h(x)dx 


( 1 ) 


The PSD, as it is defined, is the “spectrum" of the surface roughness. Its value at / is simply the power 
at that frequency. It is easy to distinguish between periodic and random rough surfaces from their PSDs. For 
periodic rough surfaces, there are some “spectral lines” in their PSDs; while these lines don t exist for a real 

random rough surface. 

Given a PSD function 2 W { , the surface roughness amplitude RMS in the frequency band of /1 - / 2 (both 
fi and / 2 are positive) can be calculated as: 
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3. CHANDRA X-RAY OPTICS 

The Chandra X-ra V optics - High Resolution Mirror Assembly (HRMA) - is an assembly of four nested Wolter 
Type-I (paraboloid and hyperboloid) grazing incidence mirrors made of Zerodur and I coated with indium (Ii )_ 
The mirror elements were polished by Hughes Danbury Optical Systems, Inc. (HDOS) m Danbury CT. The 
surface roughness was measured during the HDOS metrology measurements after the final polishing, but before 
the iridium coating. 7 Tests conducted on sample flats before and after the coating indicate that the coating 
does not change the surface roughness. 

The instruments used for the measurements were the Circularity and Inner Steticn (CIM), the 

Precision Metrology Station (PMS), and the Micro Phase Measuring Interferometer (MPMI aka WYKO). The 
CIDS was used to determine the circularity and the inner diameters. The PMS was used to measure along 
individual axial meridians. With these two instruments, HDOS essentially measured the ‘hoops and 'staves of 
each mirror barrel, and thus mapped the entire surface. The micro-roughness was sampled along meridian a 
different azimuths using the WYKO instrument at three different magnifications (* 1 . 5 , *10 x 4 U). ’ 

These metrology data were Fourier transformed and filtered. The low frequency parts of the CIDS and PMS 
data were used to form mirror surface deformation (from the designed mirror surface) maps. The high frequency 
parts of the PMS data and the WYKO data were used to estimate the surface micro-roughness. Both of them 
are parts of the HRMA model we built for the raytrace simulation of the Chandra performance. 

*The definition 2 Wi is conventional, where the subscript 1 denotes 1 -dimensional, the PSD satisfies PSD(-f) - 
pSD(f ), and typically positive frequency limits are used for most spectral integrals. The total power, <7-, is the integral 

of 2 W\ from / = 0 to oc, i.e. o~ = J 0 2 W i (f)df. 


Table 1. HRMA Mirror Sections and Their Surface Roughness 
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1.95 

2.07 
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61 


The mirror surface micro-roughness has little variation with azimuth, but tends to become worse near the 
mirror ends. We divided the data for each mirror into several axial sections which were selected so that the 
measured roughness at several places within a section were reasonably uniform; this resulted in a total of 61 
sections. We then averaged the PSD measurements within each section to provide an estimate of the PSD 
for that portion of the mirror element. Table 1 shows the resulting surface roughness in the 61 HRMA mirror 
sections. The eight mirrors are named Pl,3,4,6 (paraboloid) and HI, 3, 4, 6 (hyperboloid) due to historical reasons 
(there were 6 mirror pairs when the HRMA was designed). The number underneath each section name is the 
surface roughness amplitude RMS, <7i-iooo/mm, calculated according to Eq. (2) for / = 1 - 1000 mm -1 . Each 
mirror is 838.2 mm in length. The middle sections (M), which are the best polished and hence have the lowest 
PSDs, cover most part of the mirror surface (the number in parentheses after each M denote the percentage 
coverage). The a s for the M sections are only 2-3 A. The end sections, where the a's are relatively higher, 
cover a very small part of the mirror (< 1%), and hence contribute very little to the mirror performance. 

Figures 1 and 2 show the PSDs of the M (middle) and SC (small end) sections of PI. PI and HI were the first 
polished mirror pair and are slightly "rougher” than other pairs (see Table 1). The dash and dotted lines show 
the data from different measurements: the PMS data are in the low frequency range (/ = 0.001 - 0.3 mm -1 ); 
the WYKO data with 3 magnifications are in the higher frequency range (/ = 0.3 - 1000 mm' 1 ). The solid 
line is the combined PSD from all four frequency ranges. The SC section obviously is much rougher than the 
M section. 


4* MODEL SURFACES 

Typically, a random rough surface is only described by its PSD. Most of the methods calculate the scattering 
from the surface PSD. However, our method calculates the scattering directly from the surface high frequency 
spatial profile. Therefore, we first need to construct a model surface that is consistent with a given PSD. From 
a random rough surface profile, one can derive a unique PSD. But from a given PSD one can’t construct the 
original surface, because the phase information was lost when deriving the PSD. However, one can construct 
many model surfaces with the same roughness as the original one from the given PSD by assigning different 
random phase factors to the spectral component. 




PSD 2W, (A 2 mm) 





To construct a 1 -dimensional model surface with length L, we need to obtain N consecutive surface height 
values hi — h(xi) with a fixed interval Ax to cover the surface (i.e. NAx = L), and its surface tangent values 
h\ = h f (xi). In Appendix A, we show that h t and h\ can be computed from the surface PSD using the following 
Fourier transforms: 
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where Hj = Ny PSDi <h ) e tiPj , \ s the assigned random phase factor, and Af = l /NAx. Both h t and h[ 
are real, this requires H-j = Hf, i.e. PSD(f-j) — PSD(fj) and <p-j = —<fj. 

To construct the model surfaces of HRMA, we choose N = 2 21 and Ax = 0.0004 mm. So L = NAx = 
838.86 mm, and Af = l/A r Ax = 0.001192 mm' 1 . Figures 3 and 4 show one set of model surface sections Pl-M 
and Pl-SC, constructed using Eqs. (3) and (4) with these parameters. 


5. SCATTERING FROM MODEL SURFACES 

In this section, we calculate the scattering of plane incident waves from a surface model; most of the detailed 
derivations of the formulae can be found in Appendices B and C. 

We assume that the surfaces are sufficiently smooth so that: 1) there is no shadowing of one part of the 
surface by another; and 2) there is no reflection from one part of the surface to another, i.e. there are no multiple 
reflections by the same surface. For an incident plane wave with grazing angle Q\, the first condition requires 
that the absolute values of all the surface tangents, |/i'|, are less than 0]. The second condition requires \h[\ less 
than 0 x/2 (when h\ = -6 1 / 2 , the reflected wave is parallel to the surface). The first condition is automatically 
satisfied when the second condition is met. So the surface smoothness condition for applying this method is: 

\h,\ < — ( 5 ) 


This condition is easily satisfied for all 61 sections of the HRMA, as can be seen by comparing the tangent 
distributions in Figures 3 and 4 with the mean grazing angles of the four shells (51. 26', 41. 27', 36.43', 27.08'), 

The scattering in the transverse direction, i.e. the off-plane scattering, is smaller than that in the plane 
of incidence by approximately a factor of the grazing angle, and consequently is less than the uncertainties 
in our surface PSD. Therefore, in this paper, we ignore the off-plane scattering and limit our discussion to a 
1-dimensional surface. 

The scattering formula is given by the discrete Fourier transform of the field on the flat surface So, as shown 
in Eq. (58) in Appendix C: 
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where the scattering intensity I is a function of the scattering angle 0j+ q / p , which is the deviation from the 
specular reflection direction towards the surface; A is the wavelength; E, is the field amplitude, after the 
reflection, on the flat surface at the uniform grid X{ where the model surface was constructed. As described 
in Appendix C.3, E; is a function of the incident wave, the model surface height and tangent, and the local 
reflectivity. A is a normalization factor given by Eq. (63). Again we choose N = 2 21 to use the entire length of 
the model surface for the FFT computation. 
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Figure 3. A model surface of Chandra mirror section Pl-M, which covers 88% of the mirror PI. The top two panels 
show the surface height and the surface tangent for a 1 mm section of the model mirror. The bottom panel shows the 
surface tangent distribution of the entire surface, which can be described accurately by a Gaussian with a 60. 
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Figure 4. A model surface of Chandra mirror section Pl-SC, which is the “worst” end-section of the PI mirror. The 
top two panels show the surface height and the surface tangent for a 1 mm section of the model mirror The bottom 
panel shows the surface tangent distribution of the entire surface, which can be described accurately by a Gaussian with 

a = 138.7". 







Field Intensity vs. Scattering Angle 

1 0 8 — • 1 1 ■ ’ r 



-80 -60 -40 -20 0 20 40 

Scattering angle 6 (arcmin) 



Encirciled Energy vs. Scattering Angle 



0.001 0.010 0.100 1.000 10.000 ioo.ooaooo.ooo 

Scattering angle (arcsec) 


Scattering Function vs. Scattering Angle 



Figure 5. The scattering of 1.49 keV X-rays at 51.26' grazing incident angle from the model surface Pl-M. The top-left 
panel shows the scattering field intensity verses the scattering angle. The very sharp peak is at the specular reflection 
direction 0 = 0. The asymmetric nature of the scattering is clearly shown. The top-right panel is the same plot but 
zoomed into the core of the peak; it shows the Fraunhofer diffraction pattern due to the finite mirror length. The 
bottom-left panel shows the fractional Encircled Energy (EE) verses the scattering angle, for both sides of the specular 
direction, and also their sum. The bottom-right panel shows the scattering function £ verses the scattering angle in the 
same range as the top-right panel. 


Figures 5 and 6 show the scattering results for 1.49 keV X-rays incident upon the mirror PI at its mean 
grazing angle ( 51 . 26 '). The top two panels show the scattering field intensity verses the scattering angle. The 
sharp peak of specular reflection (top-left) and the Fraunhofer diffraction pattern (top-right) are shown as 
expected. The bottom two panels show the fractional Encircled Energies EE + , EE - , EE and the scattering 
function 5 defined as: 
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Figure 6. Scattering from model surface Pl-SC. It has much broader scattering peak than Pl-M. 
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where £ h £ s and 7 Z are the total incident and scattered energy, and the reflectivity of the rough surface as 
described in Appendix C.5. 


6. SUMMARY AND FUTURE WORK 

We have developed a method to model the wave scattering from random rough surfaces. Model surfaces with 
the same roughness as the actual surface are constructed from the actual PSD. The scattering from the model 
surfaces is calculated using the scattering formulae we derived in this paper. These scattering formulae are 
based on the general Ivirchhoff equations but without small angle approximations. This method treats the 
reflection and scattering together and provides the dependence of the reflectivity on the surface roughness 
is applicable in general and is especially useful for X-ray scattering at grazing angles. We have applied this 
method to the mirrors of the Chandra X-ray Observatory and have shown that the calculated scattering profile 
is as expected, including the Fraunhofer scattering patterns which result from the finite length of the surfaces. 

This work is still continuing. Next we will generate scattering tables, which are the tabulations of the 
scattering function S. Then we will use these scattering tables in our raytrace model to simulate the CXO 
performance and compare it with the real results of the CXO, from both on-orbit observations and its ground 
calibrations. This method should be useful for other X-ray telescope missions as well. 





APPENDIX A. CONSTRUCTION OF MODEL SURFACES 
A.l. Fourier Transform 

The Continuous Fourier Transform equations are 9 : 


Hjf) = 

[ hjx) e‘ 2,rx/ dx 

J — oc 

( forward ) 

(ii) 

hjx) = 

/•OO 

/ Hjf) e ~ t2nxf df 

J — oc 

(inverse) 

(12) 


Here if h is a function of position, x , in mm, H will be a function of spatial frequency, /, in mm 1 . 

When there are N consecutive sampled values at x — x x with the sampling interval A.r, we make the 
transform: 


x => Xi = i Ax, 

h(x) h t = h(xi), 

t N 
l = — ( — 

v 2 

-1),... 


N 
’’ 2 

(13) 

f => fj = j A/, 

HU) =, H, = SM, 

. _ ,N 
J ^ 2 

-1),.. 

.,-1.0,1,.. 

N 

’’T 

(14) 


where A xAf = l/N. We obtain the Discrete Fourier Transform equations: 


N/2 



Hj = Y. k 

»=-(N/ 2-1) 

N/2 

** = b E *>'-'** 

j=-(N/2-l) 

( forward ) 

(15) 

(inverse) 

(16) 


A. 2, Surface Height 

From Eq (1), we obtain: 


•> r L/i 

PSDjf) = j / e ,2nxf hjx)dx 


j PSDjf) L 

r L/1 

/ e t2 ” xf hjx)dx 

L J-Lf 2 


V 2 

J-L/2 


Here PSD(f) is a real continuous function of the spatial frequency /. We first need to convert Eq (17) to a 
discrete Fourier transform. Using the equations in A.l and relation L — N Ax = 1/A/, we obtain: 


Wj\ 


_ \H{fj)\ _ 1 PSDjfj) L 


A.r 


- N i 


PS DU,) IS I 


N/2 

]T h t e ^ 

i=-( A72-D 


(18) 


Therefore Hj can be expressed as the forward Fourier transform of hi as 


Hj = N 


PSDjfj) A / 


N/2 

E 

i= — (7V/2- 


hi e 




Hence the surface height, h(xi) — h tl can be expressed as the inverse Fourier transform of Hj 


(19) 


hi = 


1_ 

N 


N/2 

E H J e ~ 

/= — (N/2 — 1) 


1 N/2 f 

-7i E 4 


PSDjfj) A / 




j = -(N/ 2-1) 


( 20 ) 


where j is a random phase factor. A set of surface heights, /i t , can be generated from a set of phase factor 
Therefore for a given PSD, we can generate as many sets of surface map (of the same roughness) as we want by 
changing the random phase factor ipj. Because h %1 the surface height, has to be real, this requires H-j ~ Hj , 
i.e. PSD(f-j) = PSD(fj) and 



A. 3. Surface Tangent 

Since 


(21) 


I v ‘ E 


o -i2nx x fj 


j=-(JV/2- 1) 


j=-(JV/2-l) 


The surface tangent can be obtained by taking the derivative on both sides of Eq. (21) with respect to x<: 
N/2 N/2 

h \ = 4 E e_,2 ’ r *‘ /j = ^ E ( 

? A/ i it 


j=-(A’/2-l) 


j=-(A 72 -l) 


The surface tangent h\ also has to be real. This condition is automatically satisfied because 

-,2 nf-jH-j = = *2jr/j Hj = (~i2x fj HjY (23) 


APPENDIX B. KIRCHHOFF SOLUTION 

The wave scattering from random rough surfaces is described by the Kirchhoff solution 2 and its far-field ap- 
proximation. 

Z 



Fieure 7. Wave scattering from a random rough surface. A flat surface So with - - 0 lies in the x-y plane (y- axis not 
, ® . r c i ias sur face height : = h(x,y ), deviates from So- The z axis is normal to the x-y plane a c 

points up Incident and reflecting (or scattering) wave-vectors are shown as k, and k 2 . Incident and reflecting grazing 
angles with respect to the surface So are 9i and 0 2 . r 0 is the observation point where the scattering is to be measured. 


As shown in Figure 7, define: 

• S 0 — 2-dimensional flat surface at z = 0. 

. s — 2-dimensional rough surface, described by its surface height s = h(x,y). 

• E,e' ki r = E, e *( fc i*+*3i) — incident plane wave (in the incident plane, therefore k> — 0). 

• Eoe ,k2 r = E>t' l{krX+k * y+k ''- ::) — reflected or scattered wave from the rough surface S. 

• 0., — incident and reflecting grazing angles with respect to the surface S 0 . 



where kj and kj are the wave vectors of the incident and scattered waves, so Ei • ki = 0, E 2 ■ k 2 = 0, and 

k = ~X = l^il = = 1^2 1 = \J k'i + (24) 

A vector normal to the local surface on 5 is given by: 

dh{x,y). dh(x,y) 


n = - V(h(x,y) - z) = 


dx 


-x — 


dy 


y + z 


(25) 


The field at an observation point ro is given by the integration of contributions from the field Ef s)e^* ,I+ * 3l l 
on the surface S : 


E( r °) = ^J s J d5E(s)e ,(fclI+t35 > ^ (n ■ r) = ^ JJ dxdyE(s)e l 


(kix+kshix^y)) £ 


ikr 


(n ■ r) 


(26) 


where ds is an element, of surface area; E(s) is given by the incident wave Ei multiplied by the suitable 
reflection coefficient; the vector r goes from the point of integration (x, y, z) to the observation point (xq, yo , zo ), 
and r = |r|; n is a unit vector in the direction of n, and (n • r) ds — (n * r )dxdy. Eq. (26) is known as the 
general Ivirchhoff solution for the wave scattering. 

Next we derive the far-field approximation of this solution. When the reflecting surface is near the origin of 
the coordinate system and the observation point is far from the origin, i.e. when (x x 0 , y yo, - we 

have: 


k 2 - k x x + k y y + k z z = k ■ — 1 - x + k -- - -- -- - y + k — ~ ^ z » — (x 0 x -f y 0 y -f z 0 z) ( 27 ) 

|r| |r| |r| r 0 

r = (x 0 x)x + (yo y)y d - (zq z)z % x 0 x + y 0 y + z 0 z « j- (k x x + k y y + k z z) (28) 

r - |r| = yj (x 0 - x) 2 + (yo - y ) 2 + (^o - z ) 2 % r 0 - — x - — y - — z (29) 

n> r 0 r 0 

where r 0 = |r 0 | = yj x$ + y^ + ~o • Keep the first order of r in the phase factor and zeroth order elsewhere: 


n r 


e 


ikr 


>'o f, dhjxjj) dh(x, y) 
k [ 1 dx +ky dy k \ 

p t*(ro- ^ e ikr 0 ^tik^x + kyy+k^hix.y)) 


(30) 

(31) 


Eq. (26) becomes: 


E(r 0 ) 


ill 


dxdy E (s)e 


I (*] X + fc3*) . 


?ikr 0 


- i(k x x+k v y+k x z ) £2 


, dh{x,y) t , dh(x,y) , 
x dx + y dy ' 


ie 


ikr 0 


2nr 0 

pikra 


If 


dxdyE(s) € *i k ^ x + k ^ h ( x >y)) e -*(k*x+k v y+k x h(x t y)) 


dh(x, y) _ ( dh(x, y) , 
x 8x + v ~~dy z 


= l l 1 [ [ dxdyE(s)e , ^ kl - k * )x ~ kvy+ ( k3 ~ kz)h{x ' y )} 

lirroJJ 


dh(x, y) t dh(x, y) 
: dx v dy 


(32) 

(33) 

(34) 


This is the far-field approximation of the Kirchhoff solution for the wave scattering. 


APPENDIX C. SCATTERING FORMULA 

In this section, we derive the scattering formula from the Kirchhoff solution for the constructed model surfaces. 


z 

A 


V 


» X 


u 



Figure 8. The scattering geometry. The Hat surface So is located on the x axis. The z axis is normal to the surface So r 
The u-v axes form a coordinate system that is rotated clockwise from the x-z axes by ( 3 " *«)■ 50 tl,e J **" . ,S **5*”® 
with the specular reflection direction. An incident ray, k,, comes ,n from the left with a grazing angle 0i; had it struck 
the surface So at xi, it would have been reflected parallel to the v axis as k 20 . However it actually strikes the roug 1 
surface S at r,(x„c,), and is reflected at an angle as k 2 . The intersection of k 2 with the surface So is at x r ,. 


C.l. Integral on 1-dimensional flat surface So 

We first reduce the Kirchhoff solution to a 1-dimensional integral on flat surface S Q . Consider: 

• In plane scattering: k y = 0 

• One dimensional surface, i.e. h(x,y) only depends on .r: h(x,y) — h(x) 


Eq. (33) becomes: 


E(r 0 ) = j dxE(s) e’ (fc ' I+ * 


-+■ k%h{x)) ^ — i{kxZ~^~k*h(x}) 


dh(x) 
r dx 


■ k. 


(35) 


here we have omitted a dimensionless factor « = Y/r 0 , where Y is the transverse surface dimension; this factor 
will be absorbed later in an overall normalization factor A. 

Figure 8 shows the scattering geometry. The incident ray, k,, strikes the rough surface at r * (*«’-») and 
is reflected as k 2 , where x, is one of the N positions of the constructed model surface (see Appendix A) and 
Zi = h(xi) = hi. The reflected field at n is 




(36) 



For the integral (35), this is equivalent to have a field at (x r , , 0 ), the intersection of k 2 and x axis, on the 
surface So described by: 


E(x r ,,0) = E{ri)e ,(k ' z ' +k3hi - khi/sin(>2) 
where khi/sinO-i is the phase delay between (x,, /i,) and (x r ,,0). Let: 

E(x r J = E(x rj , 0) e ~ ,k ' x '> = E ( ri ) e '(i'iZi+ksh i -kh,/sin0 2 -k iX r i ) _ EfrOc* 0, 

So the integral (35) can be written as 


(37) 


(38) 


E(ro) 



e i{kyx-k z x-k z h(x)) 


' dh{x) 
fox 0 

OX 


(39) 


Now the integration boundary has changed from E(s) on the rough surface S to E(x) on the flat surface 5 0 , 
so h(x) = 0 and dh ^ ^ = 0. Therefore Eq. (39) becomes: 


jpikro r JL % kr 0 r 

E(r 0 ) = E(k x ,k z ) = ^ I dxE(x) e i{kl ~ kz)x (~k z ) = z — — I dxE(x) e l{k '~ kz)x (40) 

here the reflected field E(x) are calculated at non-uniformly distributed, discrete points x = x ri . The position, 
and the phase, <p t , of the field E(jr r .) are: 


X r. — 


hj 

tan do 


(pi — k\Xj + k$hi — ^ kix rt = k ( cosO\Xi — sin0\hi cosd\ 


= —khi 


sin do 

1 


sin 0o 


( Xi - -^—)) 

V tan 8 2 ) ) 


v, + — — - r—-) = = _ 2lhi 


iin 9o tan do 


sin 60 


sin- 2 ^2 
sin do 


(41) 


(42) 


where k 3 — — ksindu because, by definition, the z axis points up. 

Thus for the field E (s) of each ray ki at r*, we can use its equivalent field E(x) at x ri to do the integral 
( x rt < Xi when hi > 0 , x rx > Xi when hi < 0 ). 


C. 2 . Fourier transform with variable f 

Define a coordinate system u-v that is rotated clockwise from the x-z axes by (f — 0i ), so the v axis is aligned 
with the specular reflection direction (see Figure 8 ). Define the scattering angfe, d, as the angle of deviation 
clockwise from the v axis, i.e. d = d\ - do. Also define the variable £ = **.,"** . Therefore: 


27t£ 

0 


k cosd 1 , k x = k cos do = k cos{9\ —9), k z — k sin do = 

d d 

ki — k x — k cosd 1 — k cos{d\ — d) — —2 k sin{d\ — — ) sin - 


d\ — cos 1 




d\ — cos 1 ( cos d\ — £A) 


k sin{d\ — 9 ) 


(43) 

(44) 

(45) 


The scattering equation (40) becomes: 

E(r.) = £«(*)) = j 


ie tkr ° sin(d\ — 9) 
A 



(46) 


Thus, the scattering field E(£) can be obtained from the Fourier transform integral of the field E(ar) on the 
surface Sq. And it can be can be expressed as E(0) using Eq. (44). 



C.3. Discrete Fourier transform at Xi 

In practice, this integral is performed numerically using the Fast Fourier Transform (FFT) on N uniformly 
distributed points Xi 's where we constructed the model surface. Therefore we need to convert the field E(a r , ) 
to the field E(aq). This can be simply done by multiplying E(x r .) with two factors: 


E(.r t ) = AiBi E(x ri ) = A t B x E(x , - -A-) = A t B,E(r,)e^ 


(47) 


Where the factor ,4, is used to adjust the incident plane wave density due to the different surface height An s 
at the uniform grid * s ’s: it is calculated by intercepting all the incident rays that strike on the surface S at 
\l M 's with a coordinate that is perpendicular to the direction of incidence. Let the intercepting points be 

w t 's on the coordinate. Then: 


Ai = 


U>j+1 - Wj-1 

2 Ax si nO i 


(48) 


The factor B t is used to adjust the outgoing ray density due to the redistribution of the reflected rays from the 
non-uniform grid x ri to the uniform grid x { . For example, when the point x ti falls between the fixed grid p , 

Xi-i and x t (x z - x,_i = Ax), then 


Xi — x 

Ax 

Xi — l 


Ax 


— E(a- r; ) is added to field E(xj-i) 
E(x r ,) is added to field E(aq) 


(49) 

(50) 


This process is done for each ray until all the fields are redistributed to the uniform grid x„ 

Having obtained the field E(*<) on uniform grid, x„ we can rewrite the scattering equation (46) as the 
discrete Fourier transform (see Appendix A.l). Let: 


x => = i A.r, E(.r) => E, = E(:r,), 


s • > 0 = i 


E(0 => Ej = 


_ E(^) 


A.r 


* = -(^-1),..., -1,0,1, 
, = -(|-l),..., -1,0,1, 


N_ 

N 

' 9 


where AxA£ = 1/A r . The scattering equation (46) becomes: 

ie lkr ° sin(0\ — 0j) 


Ej - ^ 

J Ax 


/V/2 

E E < e “ 

(N/2 — 1) 


(51) 

(52) 

(53) 


where 


Ej _ - Ai Bi E(r ; ) e 1 ®' = ,4; B, Ei R{8\ + tan ‘(/i'))e 




(54) 


where R(6x + ton' 1 (h\)) is the reflection coefficient of ray i with the local grazing angle, 0, + /on" 1 (/*'), on the 
rough surface S. h\ is'the local surface tangent of the model surface. 


The scattering intensity, I. is given as a function of the scattering angle, 0, bv. 


I(0j) = /(£(0j)) = A E(fj)E*(£j) = A 


Ax sin ( 6 1 - 0j) 

A 


■y 


N/2 

E 4 e'-V 

i=-(N/ 2-1) 


(55) 


where A is a normalization factor which we will derive in section C.5. 



C.4. Scattering formula - the Fraunhofer diffraction pattern 

With the Eq. (55), it seems that we can finally obtain the profile of scattering from the rough surface. However, 
this is not quite true, because of the discrete Fourier transform. The main disadvantage of the discrete Fourier 
transform is (what else?) ‘'discrete’'. Its shortcomings are displayed perfectly in this case. Eq. (55) is correct, 
but all of the points except the central peak (9j = 0) are calculated in the valleys of the Fraunhofer diffraction 
pattern at: 


h = - 


3 A 


3 A 


N Ax sin 9\ 


L sin 9 \ 7 


j = ±1,±2,±3,... 


(56) 


where L is the surface length. In case of a perfect surface, Eq. (55) gives l(0j) = 0 except for one point at j = 0, 
and the correct diffraction pattern from the finite surface length is not obtained. To get the diffraction patterns 
at angles between 9j and 0j+ 1, we divide 9j + 1 — 9j into p equal spaces. The diffraction pattern at 9j+ q / p (q < p) 
can be calculated as: 


7(0 


j+Q/pt “ 


A £ AxsinjO i -0j+ Q /p) y 
^ { Axsin(9i - 0j+ q/p ) y 


N/2 

E 

i=-(7V/2-l) 
N/2 


_ _ 2 *i(} + <i/p) 

E i e l w 




i=-(N/ 2-1) 


(g = 0,l,2,...,p-l) (57) 


(58) 


So instead of one Fourier transform equation on E t , we need do p Fourier transform equations on E t e l?w )v /p . 
Usually, p = 8 is sufficient to calculate very nice Fraunhofer diffraction patterns. Eq. (58) is the final scattering 
formula. It maps the field on the surface, E(x), to the field intensity of scattering, 1(9). 

C.5. Normalization 

Now let’s derive the normalization factor A introduced in Eq. (55). Let e be the energy carried by each of the 
-V incident rays of the plane wave Ei. The total incident energy, total reflected energy on the surface, <f r , 
and the total scattered energy, £ tS , are: 


£i = Ne 

N/2 N/2 

£t = E l E 'l 2 = ‘ E B i \R{6i + tan ~ l (h'i))\ 1 

i=-(N/ 2-1) i= — (JV/2 — 1) 


€. = j del ( 0 ) = A j dt |E(0| 2 

Define the reflectivity of the rough surface as: 


ft = 


c 

c. r 

Tt 


N/2 


= T7 E A 2 i Bj\R(9i+tan- l (h^))f 


i=-(A'/2-l) 


Let £ r = S s . We obtain: 
A = 


e EZ / -(JV/2- 1) A i B i 1^1 1 ^ h i) 

MIe(oF 


sNK 


£ t n 


I d£ |E(£)| 2 /df|E(Ol a 


(59) 

(60) 

(61) 

(62) 


(63) 
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